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Abstract 

We derive and analyze an Asymptotic-Preserving scheme for the Euler-Maxwell 
system in the quasi-neutral limit. We prove that the linear stability condition on 
the time-step is independent of the scaled Debye length A when A — t- 0. Numerical 
validation performed on Riemann initial data and for a model Plasma Opening 
Switch device show that the AP-scheme is convergent to the Euler-Maxwell solution 
when Ax/A — )• where Ax is the spatial discretization. But, when \/Ax — )• 0, the 
AP-scheme is consistent with the quasi-neutral Euler-Maxwell system. The scheme 
is also perfectly consistent with the Gauss equation. The possibility of using large 
time and space steps leads to several orders of magnitude reductions in computer 
time and storage. 
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1 Introduction 



The goal of this paper is to derive, analyze and validate a new Asymptotic-Preserving 
(AP) scheme for the Euler-Maxwell (EM) system of plasma physics in the quasi-neutral 
limit. The Euler-Maxwell system provides a fluid description of a plasma interacting with 
an electromagnetic wave. In the one-fluid setting where the plasma ions are supposed 
immobile (sections H] to SJ, the electron fluid obeys a system of isentropic gas dynamics 
equations subjected to the Lorentz force. The electromagnetic field is a solution of the 
Maxwell equations coupled to the fluid equations through the electrical charge and cur- 
rent. In the two-fluid case (section [5]), each electron or ion species obey its own system 
of isentropic gas dynamics equations. The restriction to the isentropic case is for sim- 
plicity only: all concepts extend straightforwardly to full Euler systems including energy 
equations. 

When scaled to dimensionless variables (see section |2]), the EM system depends on 
the scaled Debye length A which is the ratio of the physical Debye length to a typical 
dimension of the system xq. The Debye length Ad is the characteristic length scale 
associated to the coupling between the particles and the electromagnetic waves and is one 
of the most important parameters in plasma physics [HI ET] . It is usually small because 
the electrostatic interaction occurs at spatial scales which are much smaller than the usual 
scales of interest. However, there are situations, for instance in boundary layers, or at 
the plasma- vacuum interface, where the electrostatic interaction scale must be taken into 
account. This means that the choice of the relevant scale Xq may depend on the location 
inside the system and that in general, the parameter A may vary by orders of magnitude 
from one part of the domain to another one. 

In the scaled EM system, A appears both in the Ampere and Gauss equations. There- 
fore, when A is very small, a quasi-neutral regime, where the local electric charge is 
everywhere close to zero, appears. Simultaneously, A <^ 1 implies that the speed of light 
is very large compared to the hydrodynamic speeds. In the limit A — )■ 0, the scaled EM 
system formally converges to a system consisting of the Faraday equation for the Magnetic 
field, of the magnetostatics Ampere equation (i.e. without the displacement current) and 
of a stationary elliptic equation for the electric field (but which is not the usual Pois- 
son equation). This system, later on referred to as the Quasi- Neutral Euler-Maxwell 
(QN-EM) system, bears analogies with the so-called Electron-MagnetoHydrodynamics 
equations (EMH) [29]. 

This paper proposes a suitable numerical scheme for both the A = 0(1) and A ^ 1 
regimes. Physically, A~^ measures the temporal and spatial frequencies of plasma oscil- 
lations and electromagnetic waves. When A ^ 1, they are very large and impose strong 
constraints on numerical discretizations. For classical explicit schemes, the time and space 
steps At, Ax must resolve these frequencies and be of order 0(A) to prevent the onset of 
numerical instabilities. For this reason, most studies are based on quasi-neutral models 
[m [25], [331 [Ml [m [121 [ini [53] . However, when A varies from one region to the other, quasi- 
neutral models lead to the wrong solution where A = 0(1). A possible way to handle such 
situations is to decompose the simulation domain and to use the full EM or the QN-EM 
models according to whether A«;lorA = 0(l) [2ni[2Il[2Sl[SIl[SIl[S21[M]. However, 
this domain decomposition approach suffers from many drawbacks. The coupling between 
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the EM and QN-EM at the interfaces is not well-defined, which questions the physical 
reliability of the any particular strategy. Additionally, the domain decomposition must 
often be updated with time, which introduces a costly mesh adaptation strategy. There- 
fore, methods which are able to handle both regimes and are free of time and space step 
constraints related to A are much more flexible, versatile and robust. This is the route 
which is followed in the present work. 

More speciflcally, we look for Asymptotic-Preserving (AP) schemes for the EM model 
with respect to the limit A — )■ 0. The AP property can be deflned as follows. Consider 
a singular perturbation problem whose solutions converge to those of a limit problem 
P° when \ (here is the EM model and P° is the QN-EM model). A scheme 
P/^ for problem P^ with time-step S and space-step h is called Asymptotic Preserving 
(or AP) if it is stable independently of the value of A when A — and if the scheme P^i^ 
obtained by letting A — )■ in P/^ with flxed [6, h) is consistent with problem P". This 
property is illustrated by the commutative diagram below: 

pA pA 



A^O 



A-5-0 



po po 



The possibility of letting A — ?■ in P/^ with flxed (5, h) implicitly assumes that the 
stability condition on [6, h) is independent of A when A — ?■ 0. This property is referred 
to as 'Asymptotic Stability'. The concept of an AP scheme has been introduced by S. 
Jin [3l] for diffusive limits of kinetic models and has been widely expanded since then 
[2lia[3|6l[71[271|30l[36l[39lll5lll8]. 

In order to achieve the AP property, a certain degree of time implicitness must be 
introduced. In section [31 we will review various implicit schemes in view of this AP 
property and show that only one of the proposed schemes does exhibit this property. 
Speciflcally, we need a fully implicit discretization of the Maxwell equations together with 
an implicit current in the Ampere equation as well as an implicit mass flux in the mass 
conservation equation. A linearized stability analysis in Fourier space shows that the 
resulting scheme is actually AP. The implicit mass-flux strategy has already been used 
for the Euler-Poisson problem [121 [ISl [HI [25] and Vlasov-Poisson problem [H [TB] and is 
also key in the large magnetic-fleld asymptotics [151 113 HE] and in the low Mach-number 
asymptotics [23]. It is a well established fact [IQ] that, in order to enforce stability of 
the hydrodynamics equations, some numerical viscosity must be added. In section [H we 
show that consistency with the Gauss equation is obtained if corresponding numerical 
viscosity terms are added to the Ampere equation. The concepts are then extended 
to the two-fluid EM model in section [5] and a numerical validation is given in section [61 
The numerical results practically demonstrate the Asymptotic-Preserving character of the 
AP-scheme. By comparison, in highly under-resolved situations (i.e. when the time and 
space steps do not resolve the fastest scales) a classical (time-explicit) scheme exhibits 
a strong instability. Implicit method have previously been proposed in the context of 
Particle-In-Cell methods for the Vlasov equation (see [9l [381 112] for the electrostatic case 
and [21 [22l HH [56] in the electromagnetic case). For hydrodynamic models, we refer to 
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[26| [T0| B9| |50] . However, few of these methods are imphcit and none has been analyzed 
in view of the AP-property. Finally, we refer to [TB] for a recent review on AP-schemes 
applied to plasma models. 



2 The one-fluid Euler- Maxwell model 
2.1 General framework 

The one-fluid Euler-Maxwell (EM) system consists of the mass and momentum balance 
equations for the electron fluid coupled to the Maxwell equations. The mass and momen- 
tum balance equations are written: 

dtn + V ■ (nu) = 0, (2.1) 
m{dt{nu) + V • (nw (g) u)) + Vp= ~en{E + ux B), (2.2) 

where n{x,t) > 0, u{x,t) G M.'^ stand for the electron density and electron velocity re- 
spectively. They depend on the space- variable x G M*^ and on the time t > 0. We denote 
by e the positive elementary charge and by m, the electron mass. The electron pressure 
p = p{n) is supposed to be a given function of n (isentropic assumption) for simplicity. 
However, the subsequent analysis would extend straightforwardly to the case where p is 
determined by an energy balance equation. The operators V and V- are respectively the 
gradient and divergence operators and u ^ u denotes the tensor product of the vector u 
with itself. We assume that the dimension d = 3 for this presentation. We have neglected 
electron- ion collisions which otherwise would introduce a friction term in (12. 2p . This term 
could be added with no change to the subsequent theory and is omitted for simplicity. 

The electric field E{x,t) eM."^ and the magnetic field B{x,t) G M.'^ are solutions of the 
Maxwell equations: 

dtB + V X E = 0, (2.3) 
c-%E - V X = -i^oj, (2.4) 
V-5 = 0, (2.5) 
V-E = eo V, (2.6) 

where eo, fio and c are the vacuum permittivity, permeability and light velocity respec- 
tively, which satisfy eo/Uoc^ = 1. Eqs (12.31) . (12. 4p and (12.61) are the Faraday, Ampere and 
Gauss equations respectively. The divergence constraints (12. 5p . (12. 6p are consequences of 
(12. 3p . (12. 4p . as soon as they are satisfied initially, which we will assume from now on. 

Finally, the electrical charge p{x,t) G M and the electrical current j{x,t) G M"^ are 
given by 

p = e{ni-n), (2.7) 
j = —enu, (2.8) 

where rii is the background ion density, which is supposed uniform and constant in time. 
Similarly, the ions are supposed steady, so that their contribution to the electrical current 
is identically zero. 
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2.2 Scaling of the one-fluid Euler-Maxwell system 



To scale this system to dimensionless units, we introduce scaling units xq, to, uq, no, Po, 
Eq, Bq, pq, jo for space, time, velocity, density, pressure, electric field, magnetic field, 
charge density and current density respectively. To reduce the number of dimensionless 
parameters, we make the following hypotheses: 

1. The spatial and temporal scales are linked by Xq = ^io^o- 

2. The velocity scale is chosen in such a way that the drift energy and thermal energy 
scales are the same: mul = Po^o ^- For convenience, we introduce a temperature 
scale To by ksTo = Po^o ^• 

3. The density scale is fixed by the uniform ion background: no = rii. 

4. The charge density scale is fixed by the number density scale by po = eriQ. 

5. The current density scale is fixed by the density and velocity scales by jo = eno^o- 

6. The electric field scale is such that the electrical and thermal (or drift) energy scales 
are the same: eEoXo = pon-o^. 

Assumptions number 1, 3, 4 and 5 are natural. Assumptions 2 and 6 guarantee that the 
inertia force, the pressure force and the electric force have the same order of magnitude. 
With these six relations, there are only three dimensionless parameters, which are: 



The first one is the ratio of the plasma velocity to the speed of light. The second one is 
the ratio of the induction electric field to the reference electric field. The third one is the 
Debye length scaled by the reference space scale. 

In this scaling, the EM system is written (by abuse of notation, we keep the same 
notations for the dimensionless variables as for the physical variables): 



We are interested in the limit A — )■ (quasineutral limit). To choose how the remaining 
parameters a and /3 scale with A, we adopt the principle of the least degeneracy, i.e. we 
choose the scaling which produces the limit system with the largest number of terms. If 
we examine fl2.13p . we notice that whatever the choice of a, we have X^a^dtE <^ a^nu. 
So, in the limit A — ?■ 0, of these two terms, only a^nu remains. The principle of least 
degeneracy thus imposes that the remaining term of fl2.13p i.e. A^/3^V x B be of the same 




(2.9) 



dtu + V ■ {nu) = 0, 

dt{nu) + V ■ {nu ®u) + Vp{n) 

P'^dtB + V X E = 0, 

X\a^dtE - P^V xB) = a^nu, 

V ■ = 0, 

A^V ■ ^ = 1 - n. 



n{E + f3^u X B), 



(2.10) 
(2.11) 
(2.12) 
(2.13) 
(2.14) 
(2.15) 



5 



order of magnitude as a'^nu, which imposes = a^. Now, the choice = 1 is the least 
degenerate one as regards eqs. fl2.1ip and fl2.12p because, either /3 <^ 1 or /3 ^ 1 will then 
lead to reduced equations with a smaller number of terms. Based on these considerations, 
we choose 

a = \, (3 = 1, (2.16) 

which leads to the final form of the scaled Euler-Maxwell system: 

+ V ■ (nV) = 0, (2.17) 

9t(nV) + V ■ (nV ® u^) + Vp(n^) = -n^{E^ + x B^), (2.18) 

dtB^ + V X = 0, (2.19) 

X'^dtE^ -V xB^ = n\^, (2.20) 

V-5^ = 0, (2.21) 

A^V ■ = 1 - n^, (2.22) 

where we have highlighted the dependence of the solution upon the parameter A. 



2.3 Quasi- neutral limit A — > 

In the limit A 0, we suppose that — ?■ n°, — )■ vP, .... Then, formally, the scaled 
EM system leads to the Quasi-Neutral Euler-Maxwell (QN-EM) system 



V ■ n° = 0, (2.23) 

+ V ■ (n° ® n°) = -(E° + n° x 5°), (2.24) 

dtB^ + V X E° = 0, (2.25) 

-V X = u\ (2.26) 

V ■ 5° = 0, (2.27) 
n° = 1, (2.28) 



The divergence free constraint on vP is a consequence of (12.261) . while the divergence 
free constraint on B^ is a consequence of (12.251) (and of the divergence free initial data). 
Finally, ra^ = 1 is no more a dynamical variable of the problem. Therefore, the core three 
equations of the QN-EM model are ([2JlD, fl2:25D . ^TM) . 

In this model, the time evolutions of vP and B^ are constrained by (12.261) . E^ is the 
Lagrange multiplier of this constraint. To resolve it and find an explicit equation for 
it suffices to take the curl of flZ^ . add it to fl^T^ and use <!d7M to cancel the 
time-derivatives. This leads to: 

V X (V X E°) + = -V ■ (m° ® M°) - M° X (2.29) 

which is a well-posed elliptic equation for E^ (provided suitable boundary conditions 
are given, such as perfectly conducting or absorbing boundary conditions ; we will treat 
the question of boundary conditions in relation to the numerical examples). In the QN- 
EM model, the hyperbolic character of the Maxwell equations is lost: E'^ adjusts to the 
variations of B^ instantaneously. 
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More precisely, the QN-EM model (E^SD-dS^H]) is equivalent to: 



(2.30) 

(E° + n°x5°), (2.31) 

(2.32) 

V- (m°®m°) X 5°, (2.33) 

(2.34) 
(2.35) 

if and only if u^\t=o and -B°|t=o are related by 

-V X B\=o = u\=o. (2.36) 

Indeed, the 'only if part of the statement has just been proved. To prove the 'if part, 
we take the curl of fl2.32p . add it to fl2.3ip and use fl2.33p to deduce that 

at(V X + n°) = 0. (2.37) 

Then, if (12.360 is satisfied, ( I2.26P is satisfied for all times. We will look for AP schemes 
which are consistent with the form fl2.30p - fl2.35p of the QN-EM model. 

If the initial conditions of the EM model do not satisfy fl2.36p . an initial layer occurs, 
during which high frequency oscillations are produced. The QN-EM model produces some 
kind of time averaging of these high frequency oscillations. The AP scheme introduces 
numerical dissipation which damps out these fast oscillations in order to approach the 
quasi-neutral dynamics. 



V -m" = 0, 

dtu^ + V ■ (u° ® n°) = - 
9t5° + V X E° = 0, 

V X (V X E°) + E° = - 

V ■ = 0, 



Remark 2.1 // we neglect the inertia of the electrons, which amounts to removing the 
drift term in the momentum equation ^2.24 ), the QN-EM model reduces to: 



dtB^ + V X E° = 0, 
^0 = -V X 5°, 
E° + M X 5° = 0, 
V ■ B° = 0, 

which is the so-called Electron-MagnetoHydrodynamics (EMH) system I2^ . Here, we do 
not make any assumption about the electron time scales, which leads to a slightly more 
complex dynamics. 

In the limit A — )■ 0, the type of the equation for the electric field changes completely, 
from a hyperbolic equation (the Ampere law fl2.20p ) to an elliptic one fl2.29p . This is 
the signature that the EM model is a singularly perturbed problem in the limit A — t- 0. 
In the process of building an AP scheme, the first step is to reformulate the problem in 
such a way that this singular perturbation character appears more explicitly. This task 
is performed in the next section. 
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2.4 Reformulation of the EM model for finite A 



In this section, we plan to find an equivalent formulation of the scaled EM model in such 
a way that the electric field equation appears as a singular perturbation of the electric 
field equation ( K29\i of the QN-EM model. With this aim, we take the curl of fl239|) . add 
it to (12. 181) . and use fl2.20p to eliminate the time derivatives of nu and B. This leads to 

X'^d'^E^ + V X (V X E^) + n^E^ = -V ■ {n\^ ® u^) - Vp{n^) - x B^. (2.38) 

In this form, it is clear that, when A — )■ and n — > 1, fl2.38p formally tends to (12. 29 p . This 
equation is a wave equation for E with wave-speed A~^. It replaces the Ampere equation 
(I2.26P in the reformulated Euler-Maxwell (REM) model: 



dtn^ + V ■ (n^u^) = 0, (2.39) 
ai(nV) + V ■ (nV ® u^) + Vp(n^) = -n^{E^ + x B^), (2.40) 
dtB^ + V X = 0, (2.41) 
X'^d'^E^ + V X (V X ^^) + n^E^ = -V ■ {n\^ ® u^) - Vp{n^) - n^u^ x B^, (2.42) 
V ■ 5^ = 0, (2.43) 
X'^V ■ E^ = (l-n^), (2.44) 



We stress the fact that this system is equivalent to the initial EM model, provided that 
E satisfies (12.261) at the initial time. This condition provides the Cauchy datum on dtE 
requested by this second order problem. 

The use of the REM model preferably to the EM model, in conjunction with an 
implicit time discretization of (12.421) . is the key for the build-up of an AP scheme for the 
EM model in the quasi- neutral limit A — t- 0. 



2.5 Linearization of the EM model 

The numerical stability analysis will use the Fourier analysis of the linearized system. In 
this section, we investigate the linearization of the EM and QN-EM models about the 
uniform stationary state n'^ = 1, u'^ = 0, E^ = 0, S'*' = 0. Expanding = 1 + en^, 
= eu^, E'^ = eE^, B^ = eB^, with e <^ 1 being the intensity of the perturbation to 
the stationary state, and retaining only the linear terms in e, we find the linearized EM 
model (in scaled units): 

dth^ + V ■ = 0, (2.45) 
dtu^ + TVn^ = -E^, (2.46) 
dtB^ + V X = 0, (2.47) 
X'^dtE^ - V X 5^ = n^, (2.48) 
V ■ 5^ = 0, (2.49) 
A^V ■ E^ = (2.50) 



8 



with T = p'{l). Introducing n^, u^, E^, B^, the partial Fourier transforms of n^, u^, E^, 
B^ with respect to x, we are led to the following system of ODE's: 

dth^ + ■ = 0, (2.51) 

dtU^ + iT^n^ = -E^, (2.52) 

dtB^ + i^x E^ = 0, (2.53) 

X'^dtE^ -i^xB^ = m\ (2.54) 

i^-B^ = 0, (2.55) 

iA^e ■ = -n^, (2.56) 

where is the Fourier dual variable to x. We denote the solution of this system by 
U^{^,t) = [n^ , , E'^ , B^) . We look for solutions of the form of a Laplace transform 
U^{iit) = e~^*f/Q(^). A simple algebra leads to the solution s = as well as to two 
non-trivial solutions: 



1. The electromagnetic mode: 



s 



it, em 



±Ul + m'^\ (2.57) 



A 

associated with the polarization E^ ± C,, 
2. The electrostatic mode: 

sles = ±{a+TX'm'^^ (2.58) 

associated with the polarization E^ \\ C,- 

In the limit A — >■ 0, both g„ and tend to oo, which corresponds to high 

frequency oscillations of the solution ?7'*'(^, t). The only mode of the QN-EM corresponds 
to s = 0. It is indeed easy to see that the linearized QN-EM model 

(9tn° + V ■ M° = 0, (2.59) 

dtu^ + rVn° = (2.60) 

dtB° + V X E° = 0, (2.61) 

-V X = (2.62) 

V ■ 5° = 0, (2.63) 

= (2.64) 

has only steady-state solutions n'^ = 0, E^ = (with adequate boundary conditions), 
while 5° is any steady-state field satisfying fl2.63p and = —V x B^. 



3 Time-semi-discretization, AP property and linearized 
stability 

We denote by 6 the time step. For any function g{x,t), we denote by g"^{x) an approx- 
imation of g{x,t"^) with t*" = mS. We present different time-semi-discretizations of the 
problem which are classified according to their level of implicitness. 
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3.1 Time-semi-discretizations of the EM system 

As mentioned in section [H we will consider different levels of time-implicitness. We recall 
that we need at least a semi-implicit discretization of the Maxwell equations otherwise the 
scheme is unconditionally unstable. As a consequence, the Lorentz force in the momentum 
equation must also be evaluated implicitly. This will be the first level of implicitness. The 
second level takes the current in the Ampere equation as well as the mass flux in the mass 
conservation equation implicitly. The third level considers a fully implicit discretization 
of the Maxwell equations, in addition to the previous levels of implicitness. 

All these schemes can be put in a unified framework by considering the following 
discretization: 



^-l^^A,m+l^A,m+l _ ^A,m^A,m^ ^ y . ® u^'^) + Vpiu^'"") = 

_^^A,m+l-a^A,m+l ^ ^A,m^A,m ^ ^X,m-j^ (^3 3) 

^2^-l^^A,m+l _ ^X,m.^ _ y ^ ^A,m+c _ ^X,m+a^X,m+a ^ ^3 ^^) 

V ■ 5^'™+' = 0, (3.5) 

.^A,m+1 _ _^A,m+l^^ (3.6) 



with a, b and c taking the values or 1. The various cases are as follows: 

1. First level of implicitness: (a, b, c) = (0, 1, 0) or (0, 0, 1): the scheme is semi-implicit 
in the Maxwell equations. The Lorentz force is implicit. The rest is explicit. This 
is the classical strategy. 

2. Second level of implicitness: (a, 6, c) = (1,1, 0) or (1, 0, 1): additionally, the current 
in the Ampere equation and the mass flux in the mass conservation equations is 
implicit. 

3. Third level of implicitness: (a, 6, c) = (1,1,1): the Maxwell equations are fully 
implicit as well as the current in the Ampere equations and the mass flux in the 
mass conservation equation. 

We note that the mass flux in the mass conservation equation and the current in the 
Ampere equation must have the same degree of implicitness in order to guarantee the 
consistency with the Gauss equation. The various schemes will be referred to by the 
value of the triple (a, b, c). For instance the (1,0, l)-scheme will refer to the scheme with 
(a, 6, c) = (1,0,1). With this level of implicitness, it is convenient to use an explicit 
evaluation of the density in the Lorentz force (13.21) . because this reduces the complexity 
of the inversion of the implicit scheme. This choice does not restrict the AP-character of 
the scheme (when applicable) nor does it change its linearized stability properties. 

We note that the first level cannot be AP. Indeed, taking the limit A — ?■ in the (0, 1, 0) 
or (0, 0, 1) schemes, we find that they do not lead to a valid recursion which allows the 
computation of the variables at time m + 1 from the knowledge of those at time m. 
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The second level could be AP. If we let A — > in the (1,1,0) scheme, we find the 
following recursion: 

^0,m+1^0,m,+l = _V X fiO'"^, (3.8) 

V ■ = 0, (3.10) 

Taking the curl of (13.71) and adding to (13. Sp . the third equation can be recast into the 
following equation for E^''^'^^: 

^0,m^O,m+l _ _y ^ ^y ^ ^0,,n^ _ ^ _ ^^0,m^O,m ^ ^0,m^ _ Vp(n°''") 

and the scheme is consistent with the QN-EM model (I2.32p - (l2.34p . It is also obviously a 
valid recursion. 

If we let A — )■ in the (1, 0, 1) scheme and we use the same computation, we find the 
following recursion: 

^0,m+1^0,m+l = _V X (3.13) 

V-5°''"+i = 0, (3.15) 

and again, the scheme is consistent with the QN-EM model (I2.32p - (l2.34p and provides a 
valid recursion formula. 

It seems that both the (1, 1, 0) and the (1, 0, 1) schemes would be good candidates AP 
schemes. However, in a forthcoming section, we will see that they are not linearly stable. 
By contrast, the (1, 1, 1) scheme will be found linearly stable. It is AP because, if we let 
A — in the (1, 1, 1) scheme, we find the following recursion: 

^-1^^0,m+l _ QO,ra^ ^ y ^ E^'^+l = 0, (3.16) 
^0,m+1^0,m+l = _V X (3.17) 

^o,m^o,m+i + V X (V X E°'"'+^) = -V • (n°'"^u°'™ (g) u^'"*) - Vp(n°'™) 

_^0,m^O,m X 5°'™, (3.18) 

V-5°''"+i = 0, (3.19) 

which is obviously consistent with the QN-EM model. It also provides a valid recursion 
for all the variables. 

3.2 Linearized stability analysis 

The goal of this section is to analyze the linearized stability properties of the previ- 
ous schemes. More precisely, we want to show that only the (1, 1, 1) scheme has the 
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Asymptotic Stability property when A — i- 0, under a suitably defined CFL condition in- 
dependent of the value of A when A — )■ 0. We will prove L^-stability uniformly with 
respect to A for the linearization of the EM model fl2.45l) - fl2.50l) . 

In general, time-semi-discretizations of hyperbolic problems are unconditionally un- 
stable. This is because the skew adjoint operator dx has the same effect as a centered 
space-differencing. For fully discrete schemes, stability is obtained at the price of adding 
numerical viscosity. To mimic the effect of this viscosity, in the present section, we will 
consider the linearized Viscous Euler-Maxwell (VEM) model, which consists of the lin- 
earized EM model fl2.45l) - fl2.50l) with additional viscosity terms (in this section, we drop 
the tildes for notational convenience): 

dtn^ + V -u^- (3An^ = 0, (3.20) 

dtu^ + TVn^ - 13 Au^ = -E^, (3.21) 

dtB^ + V X = 0, (3.22) 

X'^dtE^ -V xB^ = u^- l3Vn^, (3.23) 

V-B^ = 0, (3.24) 

A^V ■ E^ = -n^, (3.25) 

where /3 is a numerical viscosity coefficient. We keep in mind that, in the spatially 
discretized case, /3 is proportional to the mesh size h: 

13 = -fh, (3.26) 

with the constant 7 to be specified later on. To keep the consistency with the Gauss 
equation, we need to add a numerical viscosity contribution into the Ampere equation. 
The time-semi-discretization of this model leads to 

^-i^^A,m+i _ ^x,m^ ^ y . (m^'"^+«) - /3An^'™ = 0, (3.27) 
^-i(^A,m,+i _ ^x,m^ ^ TVn^'™ - /3Am^'™ = (3.28) 

^-l(^A,m+l _ ^X,m^ ^ y ^ ^A,m+6 _ ^g ^g) 

^2^-l^^A,m+l _ ^X,rn^ _ y ^ ^A,m+c _ ^X,rn+a _ f^^^\,m^ ^g gQ) 

V ■ = 0, (3.31) 

. E^'^'+i = (3.32) 

Passing to Fourier space with being the dual variable to x, we find the following recursion 
relations: 

^-i(^A,m+i _ ^x,m^ ^ iT^h^'"" + V'"* = -i^'^+i, (3.34) 

^-l(^A,m+l _ j^X,m^ ^-^^ ^X,m+b ^ ^3 35^, 

^2^-l^^A,m+l _ ^X,m^ _ ^ j^X,m+c ^ ^X,m+a _ (3.36) 

li ■ 5^'"+^ = 0, (3.37) 
^^2^.^A,m+l _^A,m+l_ ^333^ 
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All solutions of this recursion can be found as linear combinations of elementary solutions 
of the form [/^'"(^ = where e C and f/^'"* = (n^'™, m^'"", B^'"") 

Elementary algebra shows that the characteristic roots q'^{C) are the solutions of the two 
polynomial equations: 

Px{q) = - l)\q - 1 + m'^) + q'^^Q ' 1 + m?^) + q^^'^^q - 1) = 0, (3.39) 
for the electromagnetic modes and 

Qx{q) = \\q - - 1 + m'S) + q'^'-'S' + T6'\'\^W + 

+pS\\q~l + m'm\' = 0, (3.40) 

for the electrostatic ones, where we have defined 

d = b + c. 

In particular, this shows that whatever choice of the semi-implicitation of the Maxwell 
equations (either at the level of the Faraday equation or at the level of the Ampere 
equation), the linearized stability properties of the schemes are the same. 

A necessary and sufficient condition for stability is that |q''*'(OI < 1- However, 
requesting this condition for all ^ G M is too restrictive. To account for the effect of a 
spatial discretization in this analysis, we must restrict the range of admissible Fourier 
wave- vectors ^ to the interval [— f , f ]• Indeed, a space discretization of step h cannot 
represent wave- vectors of magnitude larger than |. This motivates the following definition 
of stability: 

Definition 3.1 The scheme is stable if and only if 

|gi(OI<l> such that 1^1 < ^- (3-41) 

Now, our goal is to find which of the schemes are stable under a sufficient conditions 
on 6 which is independent of A when A — (Asymptotic Stability). We prove: 

Proposition 3.2 (i) The schemes (0,1,0), (0,0,1), (1,1,0), (1,0,1) are not Asymptoti- 
cally Stable. 

(a) The scheme (1, 1, 1) is stable under the CFL condition 6 <rh where T is a constant 
independent of A and is therefore Asymptotically Stable. 

Proof: (i) Let us examine the (0,1,0) and (0,0,1) schemes first, i.e. with a = and 
d = 1. In either cases, the polynomials f l3.39p . fl3.40p can be written: 

P,iq) = X^PM + Po(g), Qx{q) = A'gi(g) + Qo{q), (3.42) 

where Pq, Pi, Qq, Qi are independent of A. More precisely, we have 

degPi = 3, degPo = 2, degQi = 2, degQo = 1, (3.43) 
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where deg refers to the degree of the polynomiaL Therefore, one of the characteristic 
roots of either equations tends to infinity and behaves hke 5^(1 + |^P)/A^ when A — 
for the electromagnetic mode and 5^/(A^(l +T5^|^p + for the electrostatic mode. 

In either cases, an instability develops when A — )■ with fixed 5. 

Let us now examine the (1, 1,0) and (1,0, 1) schemes, i.e. with a = 1 and d = 1. In 
either cases, we have 

degPi=3, degPo = 3, degQi = 2, degQc = 2, (3.44) 

Let us consider the electromagnetic mode. The two roots of Pi are 1 and 1 — 
None of them is a root of Pq as soon as > 0. In these conditions, it is easy to see 

that the roots of f l3.39p are continuous with respect to A as A — 0. Their limit q^{C) is 
therefore a solution of Po{q) = 0, which is a cubic equation with obvious root g = 0. The 
two remaining roots are easily found to be 



When 1^1 is large, the negative root becomes less than —1, which implies instability of the 
scheme. Since, when the space step h ^ 0, the maximal admissible wave-vector tends 
to infinity, there is no hope to counter-balance this instability by any restriction on the 
numerical parameters. 

(ii) For the (1, 1, 1) scheme, we have a = 1 and d = 2. For the electromagnetic mode, we 
use the same method as for the case a = 1 and d = 1, but now g = is a double root of 
Pq and the remaining root is: 

We always have q < 1 and g > if and only if 5 < (1 + |'CP)/(/5|'Cl^)- With the condition 
1^1 < 7r//i, and fl3.26p . a sufficient condition for stability in the limit 5 — > is 

1 /i^ 2 

S < -^h{l + ^)< -^h, (3.45) 

under the additional restriction h < ir which can always be assumed. For the electrostatic 
mode, a similar strategy can be developed and we notice that g = is a double root of 
Qq and no additional stability condition is required. Now, by the continuity of the roots 
with respect to A, there exists F with < F < 2/(77?^) and Ao(F) > such that under 
the condition 6 < Th, \^\ < n/h, and A < Aq, all characteristic roots q^ satisfy \q^\ < 1. 
This proves the Asymptotic stability of the scheme. ■ 



4 Spatial discretization: enforcing the Gauss law 
4.1 One-dimensional framework 

We now concentrate on the (0, 0, 1) scheme (further on referred to as the 'classical scheme') 
and the (1,1,1) scheme (the 'AP-scheme') and we investigate the spatial discretization. 
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A specific attention will be devoted to the enforcement of Gauss's law. For the sake 
of the exposition, we restrict ourselves to the one-dimensional case. In this case, all 
unknowns of the problem only depend upon a one- dimensional spatial coordinate x G M. 
The electric field has a longitudinal component E^, and a transverse component. We 
assume a rectilinear polarization, and choose an orthonormal reference frame (e^, ey,ez) 
such that Cx is in the x direction and Cy is in the transverse electric field direction. The 
magnitude of this transverse component is denoted by Ey. Finally, the magnetic field is 
aligned with and its magnitude is denoted by B^. By the divergence free condition, 
the X component of B must be uniform, and we assume that it vanishes completely. The 
velocity has components in both the x and y directions, called Ux and Uy. 
In this geometry, the dimensionless EM model is written: 



9tn^ + 9,.(nX) = 0, (4.1) 

dtin'u^) + dx{n\u>^f + p{n')) = -n\E^ + u'^B^), (4.2) 

9t(nX) + dx{n\X) = -n\E^ - u'^B^), (4.3) 

dtB^ + dxE^ = 0, (4.4) 

X^dtE^ = n\l, (4.5) 

X-'dtE^ + dxB^ = n\l (4.6) 

X^dxE^^ = 1 - n\ (4.7) 

The associated QN-EM model is obtained by taking A — ?> 0. We get: 

El + ulBl = 0, (4.8) 

dtul = -El (4.9) 

dtBl + dxEl = Q, (4.10) 

M°=0, (4.11) 

dxBl = ul, (4.12) 

n° = 1, (4.13) 



Taking the x-derivative of fl4.10p . adding to (14.91) and using fl4.12p leads to 

-dlEl + El = Q, (4.14) 



Conversely, the QN-EM model obtained by replacing (14.121) by (14.141) is equivalent to the 
original one provided that 9a;i?^|t=o = 'u^\t=o- The proof is similar to the full 3D case in 
section 12. 3[ 

The time discretization of the one-dimensional EM model is given by (omitting the 
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exponent A): 



^~l^^m.+l _ ^m.^ ^ = 0, (4.15) 

= + n"'u'^B^), (4.16) 
5~i(n'^+iM™+i - n™M;7) + 9^(n™M>;;^) = -(^"^+1-"^;^+^ - n'^<5f ), (4.17) 

r^(5™+i - 5f ) + = 0, (4.18) 

A^^"^ (E™+i -E^)= (4.19) 

^2^-1 ^^m+i _ ^ = n"^+"M;;^+^ (4.20) 

X'^d^E^^^ = l- n'^+\ (4.21) 



again, with (a, 6, c) = (0, 0, 1) (classical scheme) or (1, 1, 1) (AP scheme). 



4.2 Spatial discretization 

Now, we introduce a spatial discretization with a uniform mesh of step h and we denote 
by Cfc the cell [{k — 1/2) /i, {k + 1/2) /i] and = kh, with G Z. Like in usual first-order 
shock capturing schemes, the fluid unknowns n and u are approximated by piecewise 
constant functions within the cell Ck and represented by cell-centered values n\^, u\]^ at 
time t"^ = m5. The electric field Ex and the magnetic field Bz are approximated at the 
interfaces 0:^+1/2 = {k + l/2)h by -E'x|'^yYi/2 ^z\]t+i/2^ while Ey is approximated by 
cell-centered quantities Ey\^. The discretization of the hydrodynamic part is performed 
by means of a first order shock capturing scheme. We denote by fn\^+i/2y fuj'k'+i/2y 
fuy\'k+i/2 numerical fluxes for the mass and x and ^/-components of the momentum 
conservation equations respectively, at time t"^ = m6 and at the cell interface Xfc+i/2- 
The fully discretized scheme is written: 

6-\n\r' -n\T)+ h-\fn\T:r/2 - /niri%) = o, (4.22) 

= -nir'^'^Exir' - {nUy)\^Bz\^, (4.23) 
^-\{nuy)\r' - inUy)\T) + h-\fuS^,/, - ^,1^1/2) = 

= -n\r'-^Ey\r' + {nUx)\^Bz\T, (4.24) 

S-\Bz\^:,), - 5.1^+1/2) + h-\Ey\^:^,' - Ey\r') = 0, (4.25) 

(^^'11+1/2 - ^^-1^+1/2) = /nir+l/2' (^-26) 

x's-\Ey\r' - Ey\T) + h~\Bz\z:,), - 5.1^1/2) = (^%)ir", (4.27) 



where 



R im ^ f D I'm- I D I"* ^ Z? |m+l /' ZT" \m+l , 771 im+l \ [a r,Q\ 



The numerical hydrodynamic fluxes /n|^i/2' /«xlfc+i/2 /tiy|fc!|_i/2 computed 
using a Local Lax-Friedrichs (LLF) scheme [IQ] (also known as the Rusanov scheme [17]; 



16 



we note that this scheme enters the class of polynomial solvers of [22]: it corresponds to 
the case of a degree polynomial). In the case (a, 6, c) = (0,0, 1) (classical scheme), the 
numerical fluxes are given by: 

/nir+i/2 = \ [{nu^)\k + (^«.)ir+i + {n\t - n\t+,)] , (4.29) 
= I [{nul+pin)) \T + {nul+p{n)) 

+f^T^i/2 {inu,)\T - inu.)\T^,)] , (4.30) 

= I [inn.Uy)\T + {nu^UyW^, + f^T^y, {{nUy)\T - (n%)|- J] . (4.31) 

In the case (a, 6, c) = (1,1,1) (AP scheme) the momentum fluxes (14.301) and f l4.3ip are 
unchanged. The mass flux fn\^^y2 is given by: 

/nir+v2 = \ [('^^.oir' + inu.)\T:i + /^r+1/2 {mt - ^ir+o] • (4.32) 

Only the central discretization part of the flux is implicit, while the numerical viscosity 
term (in factor of ^^P^ explicit. The tilde is there to make a typographic 

distinction from the explicit flux fl4.29p . Indeed, using the momentum balance equation 
fl4.23p . we can relate the implicit flux fl4.32p to the explicit one fl4.29p by the following 
relation: 



r 

jn\k+l/2 ~ Jn\k+l/2 ^ 

-^(Ajr+3/2-/.jr-i/2)-^ [(r^%)ir^.ir + (^«.)ir+i^.ir+i] • (4.33) 

This flux involves an average of El^~^^ over three neighbouring mesh points which is too 
diffusive and poorly accurate. In order to reduce numerical diffusion, we replace f l4.33p 
by the following expression: 

r 

f \m+l _ f \m _ („\m _i_ „\m\ rp |m+l 
Jn\k+l/2 ~ Jn\k+l/2 2^1'=+! )^^\k+l/2 

-^(/«jr+3/2 - - + (^s)ir+i)i?.ir+i/2- (4.34) 

This implicit flux can be viewed as an order 0{6) modification of the explicit flux. This 
simple modification is crucial in making the scheme AP. 

We now specify the numerical viscosity scheme, the quantity 

is an evaluation of the local maximal wave speed at the interface Xk+i/2- It is computed 
as follows: we introduce 



\in 7-1 im.+l I _L_ r,\'^\ TP im+l I 771 |m+l 

'^Ifc+l -^2:|fe+3/2 + \P'\k+l + ^\k ) ^x\k+l/2 + J^\k 



^"-Ifc -^a;-|fc-l/2 



V 



1^+1/2 = max (z/* 1^+1/2, z/* 1^+1) and v 1^+1/2 = min (i., 1^+^/2, i^X) 



m 
x)\k 5 



where and (respectively i^*|^i/2 ^*\k+i/2) denote the largest and smallest 
characteristic speeds of the hydrodynamic systems associated to the state (n|™, [nu 
('^'"y)lr) (respectively to the state {nl'^^^/^, {nu^)\'^_^_-^^^, i^'^y)\'k+i/2)^ ^i^^ 
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and similarly for ('T-Mj^) 1^^1/2)- Then, 

/^fc+1/2 = { I ^~^\k+l/2 I ' I ^ lfc+1/2 I } • 

The time step S must satisfy the CFL condition 6 < (max^g^ ^ ensure stability 

of the hydrodynamic part of the scheme. 

An important feature of the scheme is that the current {nu)x in the x-component of 
the Ampere equation (14.261) is evaluated by using the mass flux /n 1^^1/2- level of 

the continuous problem, these two quantities are identical. Therefore, this approximation 
is consistent. However, using the mass flux rather than the current allows us to guarantee 
a perfect consistency with the Gauss equation. Indeed, taking the difference of (14.261) 
evaluated at Xk+1/2 and Xk-1/2 and using (14.221) . we easily check that: 

\'h-\Ex\^^l^, - + = A^/^^'(i?.|r+i/2 - ^.1^1/2) + (4.35) 

We deduce that the Gauss equation is exactly satisfied at any time, provided that it is 
exactly satisfied at initialization. In the ^/-component of the Ampere equation ( 14.26p . 
the current is evaluated using the usual approximation {nuy)\^''~^"' because, in a one- 
dimensional problem, the y-component of the mass flux is independent of x and does 
not enter the mass balance. In a 2 or 3-dimensional problem, one should evaluate all 
components of the current using the corresponding components of the mass flux, to ensure 
consistency with the Gauss equation. 

We now consider the sequence of updates for the two schemes separately. 



4.3 Classical scheme (a, 6, c) = (0,0,1): time update 

In this case, the time update goes as follows: first the mass conservation eq. (I4.22p is 
used to compute n*""*"^. Then the Faraday eq. (I4.25P allows us to find B^^^, immediately 
followed by the Ampere eqs. ( lOHjl . (Km to find ^xl^+i/s and Finally, with the 

momentum balance eqs. (I4.23p . (I4.24p . we find the values of (nux)]]^^'^ and {nuy)\^~^'^. 



4.4 AP-scheme (a, b, c) = (1, 1, 1): time update and AP character 

The time update follows a different sequence. We first solve for the implicit Maxwell 
equations. We begin by computing E^^^ . To this aim, inserting (I4.25P and (I4.24p into 
dMH) to eliminate 5^|^+i/2, and {nuy)\^^ respectively, we find that (I4.27P is equivalent 

to 

(A^ + 5Mt)Ey\T' - S'h-\Ey\^+,^ - 2Ey\^+' + = ^'Ey\T + S{nUy)\T 

-6h-\B,\^^,/, - 5,1™ ,/2) - 6'h-\UX+i/2 - fuS-1/2) + 6\nu.)\^B^^, (4.36) 

or, using (I4.27P again, between time steps and t™, to 

(A^ + 5Mk)Ey\T' - 5'h-\Ey\^:i - 2Ey\r' + Ey\-^,') = X'i2Ey\- - Ey\r') 

-5'h-\UX+i/2 - fuX-1/2) + 5\nu,)\^B,\^. (4.37) 
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This equation is clearly consistent with the reformulated Ampere eq. f l2.42p . Both f l4.36p 
and fl4.37l) are discrete elliptic equations for Ey\^^^ which are well-posed provided that 
suitable boundary conditions are defined. 

We now examine the computation of -E™"*"^. We insert the expression fl4.34p of the 
mass flux into fl4.26p . This yields: 

2 m m m+1 2 m 

-^(/.jr+3/2 - Ajri/2) - + 5.ir+i/2- (4.38) 

This expression provides an explicit evaluation of Ex\^]^li2- 

Once E"^^^ and E^~^^ are known, we can compute i?^"*"^ using (14.251) . then u^~^^, u^^^ 
and n""+^ using g23]), fOil) and fH:22D respectively. 

Finally, we show that the fully discrete scheme is AP. Indeed, when A — )■ 0, fl4.35p . 
gives n\^^^ = ''^IT exactly. Since we assume consistency with the Gauss equation at time 
t = 0, which, in the case A = 0, amounts to assuming that n|2 = 1 for all /c, we deduce 
that = 1 for all k and m. Then, the remaining equations yield, in the limit A — ?■ 0: 

r 7^ I m+1 / I ™ ( f \m f \m \ "IRI™ 

0-C/x|fe_|_i/2 — /n|fc+l/2 cy \Jux\k+3,/2 Jux\k-l/2) r,l"ylA: % I fc+1 J -"2 1 A;+l/2 ' 



(4.39) 



diij^lfc -d/l (iiy|fc_^i - 2i^j^|fc +Ey\^_^} 



lfc+1/2 -^2|fc-l/2y' "T "i/lfc 

-^/^-'(Ajr+1/2 - /.jri/2) + (4.40) 
'^''(^^ir+v2 - ^^ir+1/2) + /^-'(^.ir+Y - ^.ir') = 0, (4.41) 



1 

1""?/^ ~ ^y\k ) ~^ Uuy\k+l/2 ~ JUy\k-l/2) 



MT:if2 + ^^iri/2) - ^vlkB.lT, (4.42) 



(4.43) 



with the fluxes 



/nir+i/2 = 2(«-ir + «xir+i), (4.44) 
/«ilfc+i/2 = 2 [""^Ifc* ''^^Ifc+i (''^^^ir ~ ""siIa^i)] , (4.45) 

Ajr+i/2 = I + («x«.)ir+i + /^r+1/2 Kir - sir+i)] • (4.46) 

Since the pressure p{l) is now a constant, it has be removed from (14.450 . because fluxes 
are defined up to a constant in space. Inserting (14.441) into (I4.39P and the result into 
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m \ 

\k~l/2) 



f02|) . we find 

I m+1 ^ ( im r) |m\, 

"xlfc — ~^l"x|fc+l — ^"xlfe t-^x|a,-_iJ + 

im of |m iQ/ I™ / 1"^ ^ 

^ UuAk+d,/2 '^Ju^\k+l/2^ '^Ju^\k~l/2 Ju:^\k-Z/2l 

"r^l"j/lfc+l-"-2|fc+l/2 "j/lfc l-'^^lfc+l/2 "T -Dz|fc-1/2J "T "'Vlk-l-^z 

= 0{h\l + S)), 
which is consistent with (14.111) . From there, we deduce that 

/.jr+l/2 = Oih^l + 6)), fuS+l/2 = 0{h). 

It follows that fl4.43p can be written 

S-\uy\r'-Uy\T) = -Ey\r' + 0{h), 

which is consistent with (14. 9p . A similar computation, inserting (I4.42p into (14.440 and 
using (I4.39P shows that 

fn\T+l/2=Oi6h'). 

Therefore, (14. lip is such that 

~ ~2^^y\k + ^2/lfc+l) -^zlfe+1/2 + )' 

which is consistent with ( 14.80 . The consistency of ( 14.410 with ( I4.10p is obvious. Finally, 
inserting ( I4.4ip into (I4.40p leads to 

h-\B,\T:i)2 - B.\TX)2) = ^y\T + 0(5), 
which is consistent with ( 14.120 . This proves that the fully discrete scheme is AP. 



5 Two-fluid case 

5.1 Euler-Maxwell system 

The two-fluid Euler-Maxwell (EM) system consists of the mass and momentum balance 
equations for both the electron and ion fluids coupled to the Maxwell equations. The 
mass and momentum balance equations are written: 

dtUi + V • {uiUi) = 0, (5.1) 
mi{dt{niUi) + V ■ {uiUi ® Ui)) + Vpi = eni{E + x B), (5.2) 

dtUe + V ■ {UeUe) = 0, (5.3) 

me{dt{neUe) + V ■ {rieUe ® Ue)) + Vpe = -en^iE + Ue X B), (5.4) 

where the indices i and e refer to the ions and electrons respectively. The meaning of 
the variables is the same as in the one-fluid case, section (12. ip . The Maxwell equations 
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fl2.3p -( l2l6|) are unchanged but the definition of the charge and current densities is now 
given by: 



p = eijii - Tie), 



(5.5) 
(5.6) 



J = e[niUi - UeU, 



where we assume for simphcity that the ions are singly charged. 

In the scahng, the same density and velocity scales for the ions and the electrons 
are chosen. The thermal energy scale is chosen equal to the ion drift energy scale i.e. 
rriiul = ponQ^ and an additional dimensionless parameter e corresponding to the electron 
to ion mass ratio appears: 



Apart from this, we use similar scaling hypotheses as in the one-fiuid case, section 12. 2|, 
and find the dimensionless two-fiuid EM model: 



where we have introduced the dimensionless charge and current densities and j^. 

The quasineutral A — )■ limit leads to the two-fiuid QN-EM model, in which only the 
Ampere and Gauss equations fl5.13p . (15.151) are formally modified: 



We note that we keep e fixed and finite. Taking the difference of (15.161) and (15.181) and 
using (I5.23P shows that the current should be divergence free: 



dtn^ + V ■ (n^n^) = 0, 

dtinfuf) + V ■ {nfuf ® uf) + Vp{nf) = nf{E' + uf x B^), 

e'mny^) + V ■ {n'y^ ® n^)] + Vp(n^) = -n'^{E' + x B^) 
dtB^ + VxE^ = 0, 

X^dtE^ -V xB^ = -j^ := -{nfuf - n^^), 
V ■B^ = 0, 

A'V ■E^ = p^ := - nl 



(5.8) 
(5.9) 
(5.10) 
(5.11) 
(5.12) 
(5.13) 
(5.14) 
(5.15) 



d,n^ + V ■ (ny) = 0, 

dtiny,) + V ■ {ny, (g) m") + Vp(ri°) = nO(E° + m° x 5°), 
dtnl + V ■ (n°M°) = 0, 

e'mny) + V ■ {ny ® u'J] + Vp(n°) = -n°(E° + «° x B% 

dtB° + V X E° = 0, 

-VxB' = -f :=-iny-ny), 

V ■ 5° = 0, 



5.16) 
5.17) 
5.18) 
5.19) 
5.20) 
5.21) 
5.22) 
5.23) 








(5.24) 
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which is consistent with f l5.2ip . Taking the difference of fl5.17p with e fl5.19p . we find 



dtf + V ■ 0° = (n" + e-^nl)E + {nX + e~'nyi) x B, (5.25) 

with 

/ = ® M° + p(n°)Id - (n°M° ® m° + (5.26) 

the current fiux. Then, and j° both satisfy evolution equations ( fl5.20p and fl5.25p ) 
and are related by the constraint fl5.2ip . is the Lagrange multiplier of this constraint. 
To find it, we take the curl of fl5.20p . subtract it to fl5.25p and use fl5.2ip . We find 

V X (V X E°) + (n° + e-2n°)E = V ■ 0° - (n°M° + e~^n°M°) x 5, (5.27) 

which is a well-posed elliptic equation for E. An equivalent form of the two-fluid QN-EM 
model is therefore obtained by replacing the Ampere equation fl5.2ip by its reformulation 
fl5:27p provided that 

-V X B\=o = -f\t=o. (5.28) 

The proof is similar as in the one-fluid case. 

We finally note that we can reformulate the Ampere equation in the original EM model 
by using a similar manipulation. The current equation fl5.25p and its associated fiux 05.260 
have the same expression at finite A. Then, taking the curl of 05.120 . subtracting it to 
05.25P (with finite A) and using 05.2ip . we find: 

X^d^E^ + V X (V X E^) + {nf + e'^n^)^ = V ■ 0^ - {nfuf + e'^n^) x B. (5.29) 

The reformulated EM model (REM) which consists of the original EM model in which 
the Ampere equation 05.2ip is replaced by 05.29P is equivalent to the original one provided 
that the Ampere equation is satisfied at the initial time. Again, our AP-scheme for the 
two-fluid EM model will be consistent with the REM model. 

5.2 Discrete equations 

We skip the step of the time-semi-discretization as it is similar as in the one-fluid case. The 
linearized stability analysis of the two-fluid model is left to future work. We provide the 
final spatio-temporal discretization in the one-dimensional setting for reference. The one- 
dimensional equations are not recalled. They are similar to the one-fluid case, but simply 
consist in a duplicate of the mass and momentum balance equations for each species, with 
the appropriate changes in the sign of the Lorentz force. The final discretization is as 
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follows (the notations are the same as in the one- fluid case): 

ri(n,ir' - + - ^1^1%) = 0, (5.30) 

= n.ir^-'^E.ir' + {n,u,y)\TB^\T, (5.31) 
6 ^ {{niUiy)\'^~^'^ — {niUiy)\^) + h '^{fuiy\T+i/2 ~ fuiy\T-i/2) = 

= n.ir^-'^i^.ir' - {n.uaTB.\T, (5.32) 

6-\n,\r' - ne\T) + h-\uT:r/2 - = o, (5.33) 

"""((rieMex)!™'^^ — {neUex)\T) + ^ ^(/«eilfc+l/2 ~ /«ex 1^-1/2) = 

= -n,|r'-E.ir' - (r^ene,)ir^.ir, (5.34) 

= -n,\r'-^Ey\r' + (nene.)ir^.ir, (5.35) 

^"'(^.ir+V2 - ^^ir+1/2) + h-\Ey\^^,' - Ey\r') = 0, (5.36) 

^(-^2:ir+l/2 ~ ^x\T+l/2) = ~(/nJfc+l/2 ~ (5.37) 

\2r-l/77i im+l 771 |m\ , L-l/R |7n+c n im+c \ 

" [Eylk -I^y\k) + n [Bz\k+l/2- Bz\k-^i/2) - 

= -i{n.u,y)\r''-{neU.e)\r''), (5.38) 

with (a, 6, c) = (0,0,1) (classical scheme) or (a, 6, c) = (1,1,1) (AP scheme) and where 
B^l^ and ^xir^^ are given by 

In the case (a, 6, c) = (0, 0, 1) (classical scheme), the numerical fluxes are given by: 

/nJfc+l/2 = 2 ~'~ ^'^^^^^-^l^l ~'~ ^*I^V2 ('^iir ~ '^ilfc+l)] 5 (5.39) 

+^^^\k+l/2 ((^.«»)lr - ('^.«.x)ir+i)] , (5.40) 

fuiy\k+l/2 ~ 2 [i^i'^i:cUiy)\i^ + ('^'j^jx^fj?/) |/fc4-l + 

+/^.ir+i/2 ((^.«..)ir - (^.«..)ir+i)] , (5.41) 

/nelfc+l/2 = 2 [('^e^e^')!™ + ('^e^ex)lr+l + /^elfc+l/2 (""elT ~ '^elfc+l)] ' (5.42) 

fu Jk+1/2 = I [{e^Ueul^ + p{ne))\T + {e\eul, + p{ne))\T+i+ 

+f^e\T+l/2 {ineUe.)\T ' (^e^e.) l^+l) ] , (5.43) 
fuey\'k+l/2 = 2 [^'^ i(^eUexUey)\T + (^^eWex^ey) |fc+l) + 

+/^e|r+l/2 ((^en.,)|r-(r^ene,)|r+l)]- (5.44) 

The numerical viscosities are computed separately for each species with the same method 

as in section \^?2\ In the case (a, b, c) = (1, 1, 1) (AP scheme) the momentum fluxes fl5.40l) . 
f l5.4ip . fl5.43p . f l5.44p are unchanged. Using the same assumptions as in section |4^ the 
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implicit mass fluxes /„- |^j^/2' /™e lfc+1/2 (where the tildes distinguish them from the explicit 
ones) are given by: 

f \m+l _ f \m _|_ im , imN p |m+l 

^ifuiJT+3/2 ~ fu,^Xk-\l2) + ^{"^i^iylXk + -^2:|fc+l/2- (5.45) 

f |m+l _ f |m _ I im , |m\ p im+l 

Jnelfc+1/2 ~ i'ielfc+1/2 2^2 v"'elfc+l ""elfc J "^a;- 1 fc_,_i/2 

^^(/we:rlr+3/2 ~ /«e:r 1^-1/2) ~ ^ ( ('^e^ej/ ) | fc" + (""e^ey ) | T+l ) -^2|fc!|-l/2- (5.46) 

With the same computation as in section we find that this scheme satisfies the discrete 
Gauss equation exactly: 

X^A-l/' p \m+\ _ rp im+l \ _ \m+l _ \m+l\ _ 

= X'h~\E,\T^y, - 1/2) - (n,|r - n,\T). (5.47) 

The sequence of updates for the classical scheme (a, b, c) = (0, 0, 1) is a simple gener- 
alization of section H73l For the AP-scheme (a, b, c) = (1, 1, 1), we first realize that f l5.38p 
is equivalent to 

(A^ + 6\n,\^ + e-Vir))^.ir' - S'h-\Ey\l^^,' - 2E,|™+' + ^.irY) = ^'Ey\^ 
-5iin,u,y)\T - in,u,y)\T) ' 5h-\B,\t^,„ - 5.1^1/2) + 

+5'^h ^((/Miylfc+l/2 ~ ^ ^/«ealfc+l/2) ~ (/wiylfc^l/2 ~ ^ ^/"ey lfc-1/2)) + 

+5'((niMi.)|r + ^"'(neMe.)|r) ^.ir, (5.48) 
or, using f l4.27p again, between time steps t'""^ and t™, to 



(A' + nn^\l: + ^-vir))^.ir' - ^'/^-'(^.irj - si^.ir' + ^.irY^ 



= A^(2E,i--E,ir') + 

+ 5 ^ ((/«ij/ lfc + 1/2 "~ ^ /Mealfc + 1/2) ~ (/««y 1^-1/2 ~ ^ fuey\k-l/2)) + 

+5'((niM,Jir + £-2KMe.)|r) B,\k- (5.49) 

This equation is clearly consistent with the reformulated Ampere eq. ( I5.29p . Both fl5.48p 
and fl5.49p are discrete elliptic equations for Ey]^'^^ which are well-posed provided that 
suitable boundary conditions are defined. Similarly, fl5.37p is equivalent to: 

(A' + y ((^.ir+i + ^"Vir+i) + (^.ir + ^-vin)) e,\t:'/2 = 

\2 jji Im J\( f I"* / I™ 

— -C/zlfc+1/2 "Unilfc+1/2 Jne\k+l/2) 

^2/^-1 m -2 m m -2 m 

"I 2 ((^"i^lfc + 3/2 ~^ /Mea:lfc + 3/2) ~ { fuij k~l/2 ~ ^ /"ei: I fc- 1/2 ) ) 

-y(((^.«..)ir + ^-'(^e«e,)ir) + ((n,t.,,)|r+i +e-2(n,«e,)ir+i))i?.ir+i/2-(5.50) 

The remaining updates are processed in a similar way as section l4!3l The proof that this 
scheme is AP, i.e. consistent with the one-dimensional version of the QN-EM model, is 
left to the reader. 
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6 Numerical results 



In this section, we provide a numerical validation of the AP-methodology. We will consider 
two different one-dimensional test problems. The first one is a simple Riemann problem, 
where the initial condition is piecewise constant with a discontinuity at the origin. Two 
different initial conditions will be used, respectively giving rise to shock and rarefaction 
waves. The second test problem corresponds to a more realistic physical situation: it is a 
one-dimensional model for a Plasma Opening Switch (POS) device. Both tests will be run 
in the one-fluid and two-fluid cases. We will see that, while the classical scheme develops 
instabilities in under-resolved situations (when the time or space steps are bigger than the 
finest time or space scales), the AP-methodology provides a consistent approximation of 
the solution of the limit quasineutral model. We will also show that in resolved situations, 
both the classical and AP-schemes have optimal order, i.e. 1/2 for discontinuous solutions 
(in the L} norm) and 1 for smooth solutions. 

6.1 Riemann problem 

The most general initial conditions for the Riemann problem are given by: 



where n; 7^ and u^i 7^ u^r- In the two-fluid case, initial conditions like fl6.ip are 
prescribed for rig, nj, Uexi Uix- In our examples though, we will make n; = = 1 (in 
dimensionless units) and assume that the initial discontinuity applies only to the velocity 
with Uxi = —Uxr- Indeed, in this configuration and in the absence of coupling with the 
electromagnetic field, the solution is particularly simple. Since there is no analytical 
solution of the system when the coupling with the electromagnetic field is turned on, it 
is easier to qualitatively interpret the results if the solution without coupling is simple. 

Indeed, in the absence of coupling, and if Uxi > 0, i.e. if the initial velocity configura- 
tion is towards a compression of the fluid, two outgoing shock waves starting at the origin 
propagate in opposite directions at the same speed and encompass a region of higher 
density at rest (i.e. with zero velocity). If, on the other hand, Uxi < 0, i.e. if the initial 
velocity configuration is that of an expansion, two outgoing rarefaction waves starting at 
the origin propagate in opposite directions at the same speed and encompass a region of 
lower density at rest. 

When turning on the electromagnetic field, we will consider two situations. In the 
first one, the initial values of Ey, and Uy are identically zero. Then, they identically 
vanish at all times and the quantities of interest are n, Ux and Ex- In the second one, we 
suppose that the initial is non-zero and uniform. In this case, non-zero values of Ey 
and Uy are generated. 

In the one- fluid case and in the quasi- neutral limit A = 0, the solution corresponds to 
a fluid at rest (i.e. Ux = 0) with uniform density n = 1. Then, the behavior of the scheme 
in the quasi-neutral limit can be compared to this analytical solution. As A decreases, 
the numerical solution should get closer and closer to this analytical solution. 

In the forthcoming simulation, the computational domain is chosen to be [—0.1; 0.1] 
and in the two-fluid case, the electron to ion mass ratio is taken to be e"^ = 10~^. 




X < 0, 
X > 0, 




x<0, 
X > 0, 
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6.1.1 One-fluid outgoing shock waves; zero initial magnetic field 

In this test case, the initial velocities are ul = +1 and uji = —1. We first investigate 
how the schemes behave as the coupling with the electromagnetic field is turned on, i.e. 
as A is gradually decreased. Figs. [1] and |2] shows how the classical and AP schemes 
behave when A successively takes the values A = 1, A = 10~^ and A = 10~^. These 
figures display the density n and momentum as functions of space at a given time 
t = 5 X 10~^ for the classical scheme (Fig. [T]) and for the AP-scheme (Fig. [2]). We 
observe that, when A = 1, the coupling is weak and the solution is close to that of the 
Euler equations with zero Lorentz force. On the other hand, if A = 10~^, the solution 
is close to the corresponding quasi-neutral limit, i.e. n = 1 and = 0. If A = 10~^, 
the Debye length is in an intermediate regime and the solution lies in between these two 
extremes. When A is small, the boundary values of the momentum are different from 
those of the initial conditions. This is because a very fast wave has crossed the domain 
and has changed the boundary values of the momentum. This change is allowed by the 
Neumann boundary conditions which are imposed on the fiuid quantities at the domain 
boundaries. Let us now compare the magnitudes of the momentum between the final and 
initial times for A = 10~^. For the classical scheme, these magnitudes are of the same 
order of magnitude (Fig. [T]), whereas for the AP scheme the magnitude at the final time 
is very small compared to that at the initial time (Fig. |2]). Therefore, the behavior of the 
AP scheme is consistent with the quasi-neutral limit while that of the classical scheme is 
not. The AP scheme therefore ensures a correct transition from the Euler shock to the 
quasi-neutral fiuid when A decreases. 

Another way to highlight the consistency of the AP-scheme with the quasi-neutral 
limit and the corresponding inconsistency of the classical scheme is to investigate how 
the results depend on the ratio Ax/ A of the space step to the Debye length. Figs [3] and 
m display the momentum and electric field (respectively) as functions of x at the same 
final time as before. The left and right-hand pictures correspond to the classical and 
AP- schemes respectively. The value of A is kept fixed at A = lO""^ but the number 
of discretization points is decreased from = 10^ to Nx = 10"^ and finally A'^,. = 10^, 
leading to correspondingly increasing ratios Ax/A = 0.2, 2 and 20 respectively. On the 
pictures, we observe that the AP-scheme provides a neat transition from a shock wave 
solution for Ax/A < 1 to the quasi-neutral uniform solution for Ax/A ^ 1. At variance, 
the classical scheme provides large magnitude momenta or electric fields, in contradiction 
to the quasi-neutral solution. However, these solutions are not correct solutions of the 
problem with finite A either, since the wave number of the oscillations of the solutions 
have nothing to do with those obtained in the resolved situation Ax/A = 0.2. Therefore, 
in the under-resolved situation, the classical scheme is neither good for the problem with 
finite A nor for the quasi-neutral limit. 

Fig. O displays the electron momentum as a function of x at the time t = 5 x 10~^ in 
the case A = 10~^, for A^^; = 100 (left figure) and = 1000 (right figure) discretization 
points and for both the classical and AP- schemes. We see that, for this value of A, 
the momentum computed by the AP scheme is identically zero for both choices of space 
discretization Ax, while that computed by the classical scheme keeps an 0(1) magnitude. 
In these cases, Ax/A have values respectively equal to 210^ and 210^, which shows the 
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ability of the AP-scheme to handle extremely under-resolved situations. 



Riemann Shock waves, classical scheme 



Rietnann Shock waves, classical scheme 
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Figure 1: One- fluid shock wave test case with zero initial magnetic field, n (left panel) 
and nUx (right panel) as functions of x at time t = 5 x 10~^ with A = 1, A = 10~^ and 
A = 10~^ computed with the classical scheme on = 10^ space cells. 
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Figure 2: One-fluid shock wave test case with zero initial magnetic field, n (left panel) 
and nux (right panel) as functions of x at time t = 5 x 10~^ with A = 1, A = 10~^ and 
A = 10~^ computed with the AP-scheme on = 10^ space cells. 

For Debye lengths A = 1 and A = 10~^, an approximate reference solution can be 
computed on a grid such that A^'^^ < A. The grid to compute this reference solution is 
made of 10^ cells. This grid would be suitable to compute a reference solution for the case 
A = 10~^, but it cannot be done at reasonable computational cost because the Courant- 
Friedrichs-Levy condition on the Maxwell equations requires too small time steps. Indeed, 
the reference solution is computed with the classical scheme. This computation is accurate 
since all physical space and time scales are resolved by the space and time steps. The 
reference solution can be used to perform a numerical convergence study for the cases 
A = 1 and A = 10"^. We compute relative errors in the norm. For instance, the 
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Figure 3: One-fluid shock wave test case with zero initial magnetic field, nu^ as a function 
of X at time t = 5 x 10~^ with A = 10~^ for the classical scheme (left panel) and AP- 
scheme (right panel). Ax/A respectively equals 0.2, 2 and 20 for Nx = 10000, Nx = 1000 
and Nx = 100 discretization cells. 




Figure 4: One-fiuid shock wave test case with zero initial magnetic field. as a function 
of X at time t = 5 x 10~^ with A = 10~^ for the classical scheme (left panel) and AP- 
scheme (right panel). Ax/A respectively equals 0.2, 2 and 20 for A^^, = 10000, = 1000 
and Nx = 100 discretization cells. 
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Figure 5: One-fluid shock wave test case with zero initial magnetic field, nur^ as a function 
of X at time t = 5 x 10"^ with A = 10"^ for both the classical and AP schemes. Ax/A 
respectively equals 2000 and 200 for A''^. = 100 cells (left panel) and = 1000 cells (right 
panel) . 



density error is defined by: 



max 


^num 


- ^ref 


1 






1 



where Uj-d is the density of the reference solution and rinum is the density of the approximate 
solution to test. The norm is chosen because of the discontinuities involved in the 
solution of the Riemann problem. It is shown in the literature that the best convergence 
rate for the numerical approximation of discontinuous solutions of conservation laws is 
obtained in the norm and that the corresponding order is 1/2, i.e. 0{\/ Ax). Such error 
indicators are applied to both the classical and AP schemes, and for the x-components of 
the momentum and electric field. These relative errors are plotted in Fig. |6] as functions 
of Ax. We can see that the classical scheme is slightly more precise than the AP-scheme, 
and that both verify the theoretical order of convergence of 0{\/ Ax). Indeed, the slope 
of the error curve is compared to a straight line of slope 1/2 and the match is almost 
perfect. This shows that the AP-scheme is consistent with the problem with finite A in 
the resolved case. 

The studies performed on this particular test case confirm that the AP-scheme is con- 
sistent with the quasi-neutral solution in under-resolved situations, and with the problem 
with finite A in the resolved situation, as an AP-scheme should do. They also show that 
the classical scheme does not capture the correct quasi-neutral regime in under-resolved 
situations. We will now confirm these trends in the forthcoming test problems. 

6.1.2 One-fluid outgoing rarefaction waves; zero initial magnetic field 

The initial velocities in this test case are ul = —100 and ur = 100. In this low 
density region appears at the center of the simulation domain, surrounded by two outgoing 
rarefaction waves. The conclusions that can be drawn from this test-case are similar as 
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Figure 6: One-fluid shock wave test case with zero initial magnetic field. Relative 
errors on n (left panel) and nUr^ (right panel) as functions of Ax at time t = 5 x 10~^ 
with A = 1 for the classical and AP- schemes. The dashed line represents the theoretical 
error, with a slope equal to V Ax. 



for the previous test-case. Figs. [7] and ID display the density (left panel) and momentum 
(right panel) as a function of space at time t = 2x 10~^ in the cases of the classical and AP 
schemes respectively, for three values of A: A = 1, A = 10~^, A = 10~^, and for Nr^ = 10^ 
cells. In this high resolution case both schemes provide the same result. For A = 1, the 
results are close to those of a simulation of the Euler equations without coupling to the 
Lorentz force. By contrast, when A = 10~^, the density is close to a uniform one but 
some oscillations are visible near the origin and still generate a large amplitude in the 
momentum variation. However, if the ratio Ax/\ is varied from values less than unity 
to large values, we observe that the AP-scheme converges to the quasi-neutral solution. 
Fig. |9] displays momentum as a function of space in the case of the classical scheme (left 
panel) and the AP-scheme (right panel), for A = 10~^ and when the ratio Ax/A is varied 
from 0,2 (i.e. with = 10^ cells) to 2 (A^^ = 10^ cells) and finally 20 {N^ = 10^ cells). 
In the last case, the momentum computed from the AP-scheme vanishes uniformly, in 
accordance with the quasi- neutral limit, while that predicted by the classical still has 
0(1) magnitude. In the intermediate case, the magnitude of the momentum predicted by 
AP-scheme is in between that obtained in the two extreme cases. 

6.1.3 Two-fluid outgoing shock waves; zero initial magnetic fleld 

We now consider a two-fiuid model consisting of electrons and ions. By contrast to the 
one-fiuid case, where only electrons are mobile, both ion and electrons are susceptible 
to bet set into motion. In this section, we investigate the ability of the classical and 
AP- schemes to describe the setting of the ions in motion. We restrict to the case of 
the outgoing shock waves with zero initial magnetic field. The initial electron density 
and velocity are taken equal to the one-fiuid case, while the initial ion density is uniform 
equal to 1 and the initial ion velocity is uniform equal to 0. Figs. ITU] and ITT] respectively 
display the densities and the momenta as a function of space, at a given time, for the 
AP schemes. The left panels are for the electron quantities, and the right panels, for the 
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Figure 7: One-fluid rarefaction wave test case witli zero initial magnetic field, n (left 
panel) and nux (right panel) as functions of x at time t = 2 x 10~^ with A = 1, A = 10~^ 
and A = 10"'* computed with the classical scheme on 10'' discretization cells. 
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Figure 8: One-fluid rarefaction wave test case with zero initial magnetic field, n (left 
panel) and nu^ (right panel) as functions of x at time t = 2 x 10"^ with A = 1, A = 10"^ 
and A = 10~^ computed with the AP scheme on a 10^ discretization cells. 
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Figure 9: One-fluid rarefaction wave test case witli zero initial magnetic field, nux as a 
function of x at time t = 5 x 10~^ with A = 10~^ for the classical scheme (left panel) and 
the AP-scheme (right panel). Ax/A respectively equals 0.2, 2 and 20 for A^^; = 10000, 
Nx = 1000 and = 100 discretization cells. 



ion ones. Three values of A are used: A = 1, A = 10~^, A = 10~^. We can see that 
the AP-scheme provides physically meaningful results. When A = 1, the electromagnetic 
coupling between the electrons and ions is weak. The ions stay immobile with uniform 
density while the electrons exhibit outgoing shock waves as if there would be absolutely 
no coupling to the Lorentz force. By contrast, in the case A = 10~^, the electron density 
converges to a uniform density equal to one, apart from a small oscillation near the origin, 
and a comparable oscillation of the ion density (the ion density scale is magnified and 
appears larger than the electron one, but the order of magnitudes are actually similar). 
The ions are set in motion in opposite directions to the electrons as they should and the 
ratio of the ion to electron momentum scales like the mass ratio as they should (since 
the densities are almost the same). In the case A = 10~^, an intermediate situation is 
observed. In the resolved case, a convergence study can be performed with respect to a 
reference solution, computed in the same way as described in section [6. 1.11 Fig. [T2] shows 
the relative errors obtained on the electron and ion densities in the case A = 1. We 
can see that both scheme are convergent. The convergence rate of the AP-scheme seems 
a little bit slower than that of the classical scheme and the magnitude of the error a bit 
larger. However, this slightly lower precision is little price to pay for the AP-character 
which guarantees a proper behavior of the scheme in the small Debye length regime. 



6.1.4 One-fluid outgoing shock waves; non-zero magnetic field 

This test-case is similar to the one- fluid outgoing shock wave test-case of section [6.1.H but 
the magnetic field at initial time is taken non-zero. This magnetic field generates a non- 
zero ^/-component of the electric field Ey which sets the plasma into motion in this direction 
and consequently, generates a non-zero ^/-component of the velocity Uy. These components 
become larger as A is decreased. The magnitude of the dimensionless magnetic field 
at initial time is taken equal to 0.2. Such a value generates a y-component of the electron 
momentum nuy which is of the same order of magnitude as its x-component nUx when the 
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Figure 10: Two-fluid shock wave test case with zero initial magnetic field. Ue (left panel) 
and Ui (right panel) as functions of x at time t = 5 x 10~^ for A = 1, A = 10~^ and 
A = 10~^ computed with the AP-scheme on = 10^ space cells. 
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Figure 11: Two-fluid shock wave test case with zero initial magnetic fleld. 
panel) and riiUix (right panel) as functions of x at time t = 5 x 10~^ for A 
and A = 10"*^ computed with the AP-scheme on A'^; = 10^ space cells. 
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Figure 12: Two-fluid shock wave test case with zero initial magnetic field. Relative 
errors on Ue (left panel) and Ui (right panel) as functions of Aa; at time t = 5 x 10~^ for 
both the classical and AP- schemes and with A = 1. 

dimensionless Debye length A = 10""^. As in the zero magnetic field case, the quasi-neutral 
limit simply provides a uniform density equal to 1 and zero velocity in both components 
Ux = Uy = 0, while both components of the electric field vanish = Ey = and the 
magnetic field is uniform and equal to its value at time 0: Bz = Bz\t=o- As this Riemann 
problem is intended to mimic a whole space problem, we choose transparent boundary 
conditions, which in this simple ID example, coincide with homogeneous Silver-Miiller 
boundary conditions. However, transparent boundary conditions suppose that there are 
no electromagnetic sources outside the domain under consideration. In the present case, 
when the acoustic waves generated by the Riemann initial data escape the domain, they 
produce electromagnetic field sources outside the domain which are not accounted for by 
the homogeneous Silver-Miiller boundary conditions. To bypass this problem, we enlarge 
the domain to the interval [—0.2, 0.2] (i.e. twice the size of the domain of the zero magnetic 
field case) and we observe the results only on the domain [—0.1, 0.1] and for times shorter 
than the time needed for the perturbations generated by the boundary conditions to reach 
this subdomain. 

Fig. [13] displays E^ as a function of space at time t = 5 x 10~^ for the classical scheme 
(left panel) and the AP-scheme (right panel) in the case A = 10~^ and for A^^, = 10^, 
Nx = 10^ and = 10^ space cells. The results are close to those obtained in the zero- 
magnetic field case. We observe that, as Ax/A increases from 0.2 (in the case = 10^) 
to 20 (in the case = 10^), the AP-scheme correctly captures that the magnitude of 
the electron momentum gradually decreases from an 0(1) value to 0, as predicted by 
the quasi- neutral limit. By contrast, the momentum produced by the classical scheme 
remains 0(1) whatever large Ax/A becomes. 

Figs, [m and US] display Ey and as functions of x in the same conditions (left 
panel: classical scheme, right panel: AP-scheme). As Ax/A increases from 0.2 to 20, the 
approximations of Ey and Bz given by the AP-scheme tend respectively to zero and to a 
constant value equal to Bz\t=o, as predicted by the quasi-neutral limit. By contrast, the 
approximations of Ey and Bz given by the classical scheme exhibit strong oscillations with 
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increasing amplitudes as Ax/A increases. These approximations are neither the correct 
solutions for the finite A problem, nor for the limit quasi-neutral problem. 

Finally, in the case where A is not too small, a convergence study can be performed. 
Fig. [16] displays the relative errors in norm on (left panel) and Ey (right panel) 
computed with the classical and AP schemes as a function of mesh size. 
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Figure 13: One-fiuid shock wave test case with non-zero initial magnetic field. E^ as a 
function of x at time t = 5 x 10~^ with A = 10""^ for the classical scheme (left panel) 
and AP-scheme (right panel). Ax/A respectively equals 0.2, 2 and 20 for = 10000, 
Nr = 1000 and N^. = 100 discretization cells. 
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Figure 14: One-fiuid shock wave test case with non-zero initial magnetic field. Ey as a 
function of x at time t = 5 x 10~^ with A = 10~^ for the classical scheme (left panel) 
and AP-scheme (right panel). Ax/A respectively equals 0.2, 2 and 20 for = 10000, 
A'^,. = 1000 and = 100 discretization cells. 



6.2 Plasma opening switch 

Plasma Opening switches (POS) are devices used in pulsed power systems to deliver large 
currents in short times compared to the rising time of generators. A POS device consist of 
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Figure 15: One-fluid shock wave test case with non-zero initial magnetic field. as a 
function of x at time t = 5 x 10~^ with A = 10~^ for the classical scheme (left panel) 
and AP-scheme (right panel). Ax/A respectively equals 0.2, 2 and 20 for = 10000, 
= 1000 and = 100 discretization cells. 
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Figure 16: One- fluid shock wave test case with non-zero initial magnetic field. Relative 
j^i Qj-j-Qj- (left panel) and Ey (right panel) as functions of Ax at time t = 5 x 10~^ 

for both the classical and AP- schemes and with A = 1. 
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a transmission line (usually a coaxial transmission line) filled with a quasi-neutral plasma. 
The plasma short-circuits the two electrodes of the transmission line and prevents power 
to be delivered to the load. However, simultaneously, the electromagnetic wave gradually 
erodes the plasma by separating the ions and the electrons. Once a gap has been formed 
in the plasma, the electromagnetic wave can cross it and the tail of the power pulse can 
be transmitted to the load. This time-contraction enables the generation of very high 
power pulses. 

A preliminary validation of the AP-scheme can be performed on a reduced one- 
dimensional model of the 2D model such as in [23]. The computational domain extends 
over 2 x 10~^m, which is twice the length of the region filled by the plasma 10~^m. The 
plasma is located in the middle of the domain. Transparent (Silver-Miiller) boundary 
conditions are imposed at the domain boundaries, to avoid including the generator and 
the load in the simulation. Indeed, part of the incident wave is reflected back to the 
load as long as the plasma short-circuits the transmission line. It is therefore necessary 
that the boundary conditions allow these reflected waves to escape the domain. A similar 
phenomenon prevails at the other end of the transmission line in the opening phase of 
the device. The quasi-neutral plasma is at rest at initial time. The initial densities of 
the ion and electron fluids inside the plasma region are equal to 1 in dimensionless units 
and their velocities are both equal to 0. Outside the plasma region, there is vacuum, i.e. 
initial densities are 0. 

We assume a smooth transition profile for the plasma density between these two areas. 
The incident electromagnetic wave is supposed to be a Transverse Electromagnetic Mode, 
characterized by a rising time t™*^ = 10~^s and an amplitude for the electric component 
E^^^ = —1.8 X lO^V (see Fig. [T7]). At time t = the wave starts from the left side of 
the computational domain. The final simulation time is t = 2.5 x 10~^s. This time is 
long enough to allow for observation of the wave impact on the plasma and the resulting 
plasma motion. Unfortunately, the simple one- dimensional setting does not allow for 
the observation of the POS opening, as this phenomenon is related to plasma motion in 
transverse direction to the transmission line, which is not accounted for here. However 
the results from the AP-scheme shown below are consistent with the expected physical 
phenomena. 
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Figure 17: Incident electromagnetic wave profile as a function of time. 

Two subsets of test-cases are performed. First, a low density POS is considered, with 
an initial density of lO^^m"^. Second, a higher density POS with an initial density of 



37 



lO^^m"^ is simulated. For both the low and high density POS the temperatures of the 
ion and electron fluids is approximately 5 eV i.e. 58 x lO'^K, and a carbon plasma (C+) 
is considered. The low density POS allows for a fast penetration of the electromagnetic 
wave in the plasma, whereas the high density POS acts like a barrier reflecting the wave 
which has more difficulties to cross it. For both the low and high density POS, the one- 
and two-fluid models will be used. 

6.2.1 how density POS; one-fluid model 

In this case, the order of magnitude of the scaled Debye length in the plasma is A = 10~^. 
Then, a grid such that lOAx < A is made of 10^ cells. This grid is fine enough to resolve 
the small space and time scales (i.e. the Debye length and electron plasma period). 
These conditions ensure that the classical scheme is stable and accurate enough (given 
the computational time constraints) to build a reference solution. We denote Ax"^^^ and 
Af^^ the space and time steps used for these computations. The electron plasma period 
is Tp = 10~^°s . However, the most severe time constraint in this problem arises from the 
CFL condition for the Maxwell equations due to the explicitness of the classical scheme. 
The reference time step suitable in these conditions is At^^^ = 10~^^s. Then, 10^ time 
steps are needed to to obtain results at time t = 2.5 x 10~^s. These reference results are 
used to check the accuracy of both the classical and reformulated scheme. 

First a convergence test is realized by comparison with the reference solution. The 
numerical errors for n and Ex are recorded for the classical and AP- schemes with meshes 
consisting of 100, 500, 1000, 5000 and 10000 cells. Figure [18] compares the relative 
error as a function of Ax on n and between the two schemes. Both show exactly the 
same error. Moreover, the slope of the error confirms that in the case of smooth solution 
both numerical schemes are first order in space. Indeed, both curves are very close to 
the theoretical error plot (dashed line with a Ax slope). In this context the classical 
scheme ensures stable computations even if the space step is much larger than the Debye 
length. Its time-step however must is bounded by the CFL condition for the Maxwell 
equations. The level of time-implicitness in the AP-scheme ensures stability regardless of 
the time-step as long as it satisfies the CFL condition of the hydrodynamic equations. 
Since both the fluid and acoustic velocities are much smaller than the speed of light, this 
provides an enormous gain in the allowed value of the time-step. 

In the following simulations, the classical scheme is used with two parameter choices. 
The first choice allows for the computation of the reference solution, as explained above. 
The second choice is space under-resolved but time-resolved. It uses a larger mesh size 
than the Debye length, namely Ax = lOA but a time step which resolves the CFL condition 
of the Maxwell equations, the fastest time-scale in these conditions as mentioned above. 
We will refer to this situation as 'under- resolved classical scheme'. The AP-scheme will 
be run in a both time and space under-resolved situation. The mesh size will be the same 
as for the under-resolved classical scheme but the time-step will be hundred times the 
time-step of the under-resolved classical scheme. 

Fig. [19] displays Ey (left panel) and (right panel) as functions of x at time t = 2.5 
ns, for both the reference, under-resolved classical and under-resolved AP- schemes. Fig. 
[20] displays E^ and nux in a similar fashion, on Fig. Fig. [191 ^6 notice that the plasma 
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prevents the transmission of the wave, as the values of Ey and at the right end of the 
plasma are almost zero. The numerical diffusion induced by the larger time-steps used 
for the under-resolved AP-scheme is noticeable, but still acceptable given the large gain 
in computational efficiency: the computing time is reduced by a factor 100. 




Ax Ax 



Figure 18: Low density POS; one-fluid model. relative error on n (left panel) and E. 
(right panel) as function of As at time t = 2.5ns. 
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Figure 19: Low density POS; one-fluid model. Ey (left panel) and (right panel) as 
functions of x at time t = 2.5ns with the reference, under-resolved classical and under- 
resolved AP- schemes. 

6.2.2 Low density POS; two-fluid model 

Both the classical and reformulated scheme behave in a similar way in the case of two- 
fluid simulations, and the same conclusions hold for the numerical convergence study. For 
instance. Fig. |21] displays UeUey (left panel) and UiUiy (right panel) in the same conditions 
as discussed for the one-fluid model. We can see that the electromagnetic wave sets 
electrons and ions into motion in the y direction in opposite directions. With a two or 
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Figure 20: Low density POS; one-fluid model. (left panel) and nur^ (right panel) as 
functions of x at time t = 2.5ns with the reference, under-resolved classical and under- 
resolved AP- schemes. 



three dimensional model where the extension in the y direction would be bounded by 
the transmission line electrodes, this would induce a segregation of the electrons and ions 
on the different sides of the transmission line. This phenomenon induces the aperture of 
the POS. In the one-dimensional situation, the densities are supposed uniform in the y 
direction and this phenomenon cannot be seen. 




Figure 21: Low density POS; two-fluid model. UeUgy (left panel) and riiUiy (right panel) 
as functions of x at time t = 2.5ns with the reference, under-resolved classical and under- 
resolved AP- schemes. 



6.2.3 High density POS; one-fluid model 

In the case of the high density POS, the Debye length and electron plasma period in the 
plasma are one order of magnitude smaller. In this situation the wave cannot penetrate 
the plasma as fast as in the low density test case. The scaled Debye length is now of the 
order of 10~^. Then, a grid with a space step such that lOAx < A (which was the ratio 
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used for the convergence study in the low density test-case) is made of 10^ cells. The 
computational cost induced by such a fine mesh is prohibitive. For this reason, we cannot 
present any convergence study in this test-case. However the fine grid used previously for 
reference is such that Ax < A and still can be used to generate a reference solution, to 
which the solution of the under-resolved classical and under-resolved AP- schemes will be 
compared. 

Fig. 122] displays (left panel) and (right panel) as functions of x at time t = 2.5 
ns, for the reference, under-resolved classical and under-resolved AP- schemes. This figure 
shows the plasma acting like a barrier on the magnetic field . In such a high density case, 
plasma waves appear at the right end of the plasma region, where an electron beam leaks 
outside the plasma. The typical wave-length of these plasma waves is 0(A), i.e. 10~^. 
Therefore, the fine grid with mesh size Ax ~ A can resolve this scale and the reference 
solution is thus able to describe these waves in a satisfactory way. By contrast, the 
coarse grid does not resolve these waves. Therefore, the under-resolved classical scheme is 
subject to instabilities generated by the impossibility of correctly describing these waves. 
The under-resolved AP-scheme does not attempt to resolve these waves, but provides 
the correct average of the oscillation and does not suffer from any instability. We notice 
the slightly larger numerical diffusion of the under-resolved AP scheme, which is the 
counterpart of the increased time-step. Still, the use of a coarse mesh size combined with 
large time-steps allows for a large reduction of the computational cost : the CPU times 
needed to compute the reference, under-resolved classical and under-resolved AP- schemes 
results are respectively 2 x 10^ s , 20 s and < 1 s. 

High density plasma opening switch High density plasma opening switch 

2 1 I I I I 1 8 1 ' ' I 

— Reference 
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Figure 22: Large density POS; one-fluid model. (left panel) and (right panel) as 
functions of x at time t = 2.5ns with the reference, under-resolved classical and under- 
resolved AP- schemes. 



7 Conclusion 

In this paper, we proposed and analyzed an Asymptotic-Preserving scheme for the Euler- 
Maxwell system in the quasi- neutral limit. The scheme is exposed in detail for a one-fluid 
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plasma model where the ions are immobile and form a fixed neutralizing background. It 
is then extended to a two-fiuid model where both ions and electrons are mobile. The 
analysis involves a proof of its 'Asymptotic-Preserving' character and that its linear sta- 
bility condition is independent of the scaled Debye parameter when the latter tends to 
zero. The numerical simulations involve comparisons between the AP-scheme to a 'clas- 
sical' scheme in the one- and two-fiuid configurations, for two different one-dimensional 
test-cases: the Riemann problem and the Plasma Opening Switch device. The numerical 
convergence study shows that both the classical and AP-scheme are convergent to the 
Euler-Maxwell solution with resolved time and space discretizations. On the other hand, 
with under-resolved time and space discretizations, the AP scheme is consistent with the 
quasi-neutral Euler-Maxwell system. Additionally, the proposed spatial discretization al- 
lows for a perfect consistency with the Gauss equation. By contrast, in under-resolved 
situations, the classical scheme leads to spurious large amplitude oscillations and insta- 
bilities. The possibility of using large time and space discretization parameters with the 
AP-scheme leads to several orders of magnitude reductions in computer time and storage. 
Future work will pursue the validation of the methodology to multi-dimensional cases and 
extend it to plasma kinetic models such as the Vlasov or Fokker-Planck-Landau equations. 
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